function [ outEFG ] = wrapperEFG( xNodes, xSamples, optEFG)
%WRAPPEREFG Returns a wrapper of the EFG aproximants and associated 
%   gradients.
%
%   INPUT:
%   xNodes      :   Set of nodes.
%   xSamples    :   Set of samples points.
%   optEFG      :   Wrapper structure that contains rho, poliBaseType 
%                   and order fields.  
%
%   OUTPUT
%   outEFG.p_samp   :   Shape functions at each sample point.
%   outEFG.dp_samp  :   Gradient of the shape functions at each sample point.


rho             = optEFG.rho;
poliBaseType    = optEFG.poliBaseType;
order           = optEFG.order;   
neightbors      = optEFG.s_near;

totalNodes      = size(xNodes,1);
totalSamples    = size(xSamples,1);
dimension       = size(xSamples,2);
baseDimension   = getPolibaseDimension(order);

samplesShapeFunction     = cell(1, totalSamples);
gradSamplesShapeFunction = cell(1, totalSamples);

for indexSample = 1 : totalSamples

    closestNodes = xNodes(neightbors{indexSample},:);
    totalClosestNodes = size(closestNodes,1);
    gradAlpha = zeros(baseDimension, dimension);
    gradN     = zeros(totalClosestNodes, dimension);
    weights   = zeros(totalClosestNodes, 1);    
    
    %   Inner loop initialization    
%     gradAlpha = zeros(baseDimension, dimension);
%     gradN     = zeros(totalNodes, dimension);
%     weights   = zeros(totalNodes, 1);    
    
%     for indexNode = 1 : totalNodes
    for indexNode = 1 : totalClosestNodes
%         distance = norm(xSamples(indexSample, : ) - xNodes(indexNode, : ));
        distance = norm(xSamples(indexSample, : ) - closestNodes(indexNode, : ));
        weights(indexNode) = gaussianWeight(distance , rho);
    end
%     B = assemblyBMatrix( xSamples(indexSample, : ), xNodes, order , weights , poliBaseType );
    B = assemblyBMatrix( xSamples(indexSample, : ), closestNodes, order , weights , poliBaseType );
%     A = assemblyMomentMatrix( xSamples(indexSample, : ), xNodes, weights , poliBaseType, order );
    A = assemblyMomentMatrix( xSamples(indexSample, : ), closestNodes, weights , poliBaseType, order );
    p = poliBase(xSamples(indexSample, : ), order, poliBaseType);
    [alpha L U]= solveAlpha (A,p);
    N = alpha * B';
    
%     ssF_aux = N(neightbors{indexSample})';
%     samplesShapeFunction(indexSample) = {ssF_aux};
    samplesShapeFunction(indexSample) = {N'};

    gradP       = gradPoliBase (xSamples(indexSample, : ), order, poliBaseType);
%     gradWeights = gradGaussianWeight (xNodes,xSamples(indexSample, : ), rho);
    gradWeights = gradGaussianWeight (closestNodes,xSamples(indexSample, : ), rho);
%     gradA       = assemblyGradMomentMatrix (xSamples(indexSample, : ), xNodes,...
%                   gradWeights, poliBaseType, order);
    gradA       = assemblyGradMomentMatrix (xSamples(indexSample, : ), closestNodes,...
                  gradWeights, poliBaseType, order);
    loadVectorGradAlpha = assemblyLoadVectorGradAlpha ...
                        (gradP(:,1), gradP (:,2),...
                        gradA(:,:,1), gradA(:,:,2), alpha );

    for d = 1 : dimension
        gradAlpha(:,d) = solveGradAlpha(L, U,loadVectorGradAlpha(:,d)');
    end

%     gradB = assemblyGradBMatrix(xSamples(indexSample, : ), xNodes, order,...
%          gradWeights, poliBaseType);
    gradB = assemblyGradBMatrix(xSamples(indexSample, : ), closestNodes, order,...
        gradWeights, poliBaseType);
    
    for d = 1 : dimension
        gradN(:,d) = gradAlpha(:,d)' * B' +  alpha * gradB(:,:,d)';
    end
    
%     gssF_aux = gradN(neightbors{indexSample},:);
%     gradSamplesShapeFunction(indexSample) = {gssF_aux};
    gradSamplesShapeFunction(indexSample) = {gradN};
end


outEFG.p_samp = samplesShapeFunction;
outEFG.dp_samp = gradSamplesShapeFunction;

end

